from BaseFunction import *


# 使用示例
if __name__ == "__main__":
    # 读取数据
    df1 = read_satellite_data("GnssData/101336")
    df2 = read_satellite_data("GnssData/101337")
    df3 = read_satellite_data("GnssData/101338")
    df = pd.concat([df1, df2, df3], ignore_index=True)
    satellite_data = df[df["SAT"] == 'G03']

    # 设置参数
    params = {
        'ele_s': 4, 'ele_d': 9,
        'azi_s': -60, 'azi_d': 105,
        'Ant_h': 7.79,
        'data_num': 25,
        'L_r': 0.19,  # GPS L1波长
        'zqb_s': 3.0,
        'peak_s': 0.3,
        'h_s': 4, 'h_d': 8,
        'method': 1,
        'p1': 8,
        'p2': 1,
        'p3': 8
    }

    # 运行反演
    results = gnss_height_inversion(df, params)

    plt.figure(figsize=(10, 6))
    for i in range(len(results['rightf'])):
        y1 = results['frp2'][i]
        detrend = results['frp6'][i]
        plt.plot(y1, detrend, '-', color=np.random.rand(3, ))

    plt.xlabel('sin(Elevation angle)')
    plt.ylabel('Detrended SNR (V)')
    plt.title('GNSS-R Reflection Signals')
    plt.grid(True)
    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
    plt.tight_layout()
    plt.show()

    plt.figure(figsize=(10, 6))

    for i in range(len(results['rightf'])):
        frp = results['frp1'][i]
        # 绘制频率变化曲线（带随机颜色）
        color = np.random.rand(3)
        plt.plot(frp[:, 0], frp[:, 1],
                 '-',
                 linewidth=1,
                 color=color,  # 生成RGB随机颜色
                 label='频率变化' if i == 0 else "")

        # 绘制峰值功率红点
        plt.plot(results['rightf'][i], results['mpeak'][i],
                 marker='o', fillstyle='none',color=color,
                 label='峰值功率' if i == 0 else "")

        # 统一设置坐标轴标签和图例
    plt.xlim(0, 100)
    plt.ylim(0, 20)
    plt.xlabel('频率')
    plt.ylabel('振幅/V')
    plt.grid(True)
    plt.legend(loc='upper left',
               handles=[plt.Line2D([0], [0], color='gray', lw=1),
                        plt.Line2D([0], [0], marker='o', fillstyle='none', color='red', linestyle='')])

    plt.tight_layout()
    plt.show()

    rel_data = read_excel_range('水位.xlsx', 'A3:B3096')
    time = [t[0] for t in rel_data]
    high = [h[1] for h in rel_data]
    time, high = interpolate_time_series(time, high)
    time1 = [t[0].replace(second=0) for t in results['bbb1']]
    high1 = [h[2] for h in results['bbb1']]

    mapping = dict(zip(time, high))
    high2 = [mapping[x] for x in time1]
    dh = [abs(h1 - h2) for h1, h2 in zip(high1, high2)]

    plt.figure(figsize=(10, 6))
    plt.plot(time, high, '-', color=np.random.rand(3, ))
    plt.scatter(time1, high1, color=np.random.rand(3, ), s=10)
    plt.xlabel('时间')
    plt.ylabel('潮位')
    plt.title('海面高度反演结果图')
    plt.grid(True)
    plt.show()

    plt.figure(figsize=(10, 6))
    plt.scatter(time1, dh, color=np.random.rand(3, ), s=10)
    plt.xlabel('时间')
    plt.ylabel('绝对误差')
    plt.title('残差序列图')
    plt.grid(True)
    plt.show()

    high1 = np.asarray(high1).flatten()
    high2 = np.asarray(high2).flatten()
    mae = mean_absolute_error(high2, high1)
    print("MAE:", mae)

    # 计算 RMSE（均方根误差）
    rmse = np.sqrt(mean_squared_error(high2, high1))
    print("RMSE:", rmse)

    # 计算相关系数 R（Pearson 相关系数）
    r, _ = pearsonr(high2, high1)  # 返回 (r, p-value)
    print("R:", r)
